namespace NBody.Analysis
{
    using System;
    using DataFile;
    
    public class AngularMomentumBin : Bin
    {
        public override void Analyse(NBodySystem s)
        {
            switch (spacing)
            {
                case Spacing.Normal:
                    bin_normal(s);
                    break;
                case Spacing.Logarithmic:
                    bin_log(s);
                    break;
                case Spacing.EqualMass:
                    bin_equal_mass(s);
                    break;
            }
        }
        
        private void bin_equal_mass(NBodySystem s)
        {
            Console.Error.WriteLine("Binning angular momentum in radial shells with equal mass spacing");
            
            int num_per_bin = s.NumParts / n + 1;
                        
            s.Sort();
            
            data = new Table(n, 8, "Radius, Mass (total enclosed), Jx, Jy, JZ, J(rms), |J|, |Jtot|");
            data.Clear();
            
            double m = 0.0, m_bin = 0.0;
            double jx = 0.0, jy = 0.0, jz = 0.0;
            double j2 = 0.0;
            double jxt = 0.0, jyt = 0.0, jzt = 0.0;
            int c = 0;
            for (int i = 0; i < s.NumParts; i++)
            {
                m += s[i].Mass;
                m_bin += s[i].Mass;
                jx += s[i].Mass * s[i].Jx;
                jy += s[i].Mass * s[i].Jy;
                jz += s[i].Mass * s[i].Jz;
                j2 += s[i].Mass * Math.Sqrt(s[i].J2);
                
                jxt += s[i].Mass * s[i].Jx;
                jyt += s[i].Mass * s[i].Jy;
                jzt += s[i].Mass * s[i].Jz;
                
                if (i % num_per_bin == 0)
                {
                    //Console.Error.WriteLine(c);
                    data[c, 0] = s[i].Radius;
                    data[c, 1] = m;
                    data[c, 2] = jx;
                    data[c, 3] = jy;
                    data[c, 4] = jz;
                    data[c, 5] = j2 / m_bin;
                    data[c, 6] = Math.Sqrt(jx * jx + jy * jy + jz * jz);
                    data[c, 7] = Math.Sqrt(jxt * jxt + jyt * jyt + jzt * jzt);
                    c++;
                    jx = jy = jz = j2 = m_bin = 0.0;
                }
            }
        }
           
        private void bin_log(NBodySystem s)
        {
            Console.Error.WriteLine("Binning angular momentum in radial shells with logarithmic spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
            double deltaR = Math.Log10(rmax / rmin) / (double)n;
        
            int[] num = new int[n];
            double[] jxbin = new double[n];
            double[] jybin = new double[n];
            double[] jzbin = new double[n];
            double[] j2bin = new double[n];
        
            data = new Table(n, 8, "Radius, Mass (total enclosed), Jx, Jy, JZ, J^2(rms), |J|, M(r)r");
            data.Clear();
            for (int i = 0; i < n; i++)
            {
                jxbin[i] = jybin[i] = jzbin[i] = 0.0;
                num[i] = 0;
            }
        
            for (int i = 0; i < s.NumParts; i++)
            {
                int pos = (int)((Math.Log10(s[i].Radius) - Math.Log10(rmin)) / deltaR);
                jxbin[pos] += s[i].Mass * s[i].Jx;
                jybin[pos] += s[i].Mass * s[i].Jy;
                jzbin[pos] += s[i].Mass * s[i].Jz;
                j2bin[pos] += s[i].Mass * Math.Sqrt(s[i].J2);
                num[pos]++;
            }
        
            double m_enc = 0.0;
            for (int i = 0; i < n; i++)
            {
                m_enc += (double)num[i] * s[0].Mass;
                data[i, 1] = m_enc;
                data[i, 0] = Math.Pow(10.0, Math.Log10(rmin) + 0.5 * deltaR + deltaR * i);
                data[i, 2] = jxbin[i];
                data[i, 3] = jybin[i];
                data[i, 4] = jzbin[i];
                data[i, 5] = j2bin[i] / ((double)num[i] * s[0].Mass);
                data[i, 6] = Math.Sqrt(jxbin[i]*jxbin[i] + jybin[i]*jybin[i] + jzbin[i] * jzbin[i]);
                data[i, 7] = 1.5 * Math.Sqrt(m_enc * data[i, 0]);
            }
        }
        
        private void bin_normal(NBodySystem s)
        {
            Console.Error.WriteLine("Binning angular momentum in radial shells with normal spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
            double deltaR = (rmax - rmin) / (double)n;
        
            int[] num = new int[n];
            double[] jxbin = new double[n];
            double[] jybin = new double[n];
            double[] jzbin = new double[n];
            double[] j2bin = new double[n];
        
            data = new Table(n, 7, "Radius, Mass (total enclosed), Jx, Jy, JZ, J(rms), |J|");
            data.Clear();
            for (int i = 0; i < n; i++)
            {
                jxbin[i] = jybin[i] = jzbin[i] = 0.0;
                num[i] = 0;
            }
        
            for (int i = 0; i < s.NumParts; i++)
            {
                int pos = (int)((s[i].Radius - rmin) / deltaR);
                jxbin[pos] += s[i].Mass * s[i].Jx;
                jybin[pos] += s[i].Mass * s[i].Jy;
                jzbin[pos] += s[i].Mass * s[i].Jz;
                j2bin[pos] += s[i].Mass * Math.Sqrt(s[i].J2);
                num[pos]++;
            }
        
            double m_enc = 0.0;
            for (int i = 0; i < n; i++)
            {
                m_enc += (double)num[i] * s[0].Mass;
                data[i, 1] = m_enc;
                data[i, 0] = rmin + 0.5 * deltaR + deltaR * i;
                data[i, 2] = jxbin[i];
                data[i, 3] = jybin[i];
                data[i, 4] = jzbin[i];
                data[i, 5] = j2bin[i] / ((double)num[i] * s[0].Mass);
                data[i, 6] = Math.Sqrt(jxbin[i]*jxbin[i] + jybin[i]*jybin[i] + jzbin[i] * jzbin[i]);
            }
        }
    }       
}
